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by 

Peter  J.  Torvik 
ABSTRACT 

A  Method  for  determining  the  two  dimensional  transient  teaperature 
distribution  and  progression  of  melting  in  disks  subjected  to  an  applied 
flux  aver  one  face  is  given.  The  calculations  are  performed  by  dividing 
the  solid  into  a  nuaber  of  finite  elements  and  performing  a  heat  balance 
over  finite  time  increments.  The  method  was  found  to  give  good  agreement 
with  known  solutions  for  two  dimensional  heat  conduction  problems  and  one 
dimensional  melting  problems.  Two  dimensional  melting  problems  in  alumina, 
titanium,  stainless  steel  and  magnesium  are  also  considered  as  examples 
of  the  method.  The  results  include  a  demonstration  that  the  tine  re¬ 
quired  to  melt  through  approaches  the  time  predicted  by  a  ono-dimensional 
(axial)  heat  balance  if  the  power  per  unit  thickness  is  sufficiently 
large. 
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I.  INTRODUCTION 

If  the  radiation  from  &  high  intensity  laser  is  directed  onto  the 
surface  of  a  solid,  significant  temperature  rises  will  occur,  and  will 
eventually,  in  the  case  of  a  continuous  incident  flux,  lead  to  melting 
and  possibly  vaporization  of  the  material  at  the  surface.  As  tine  pro¬ 
gresses,  the  region  of  melting  may  be  expected  to  increase  both  radially 
from  the  beam  center,  and  into  the  solid.  The  heating  of  the  material 
will  also  lead  to  thermal  expansion  which  will,  in  turn,  give  rise  to 
thermal  stresses  which  may  be  of  sufficient  magnitude  to  cause  fracture. 
This  phenomenon  is  most  likely  to  occur  in  materials  of  limited  ductility, 
as  ceramics;  however,  in  a  composite  material,  or  in  a  composite  struc¬ 
ture,  the  differing  coefficients  of  thermal  expansion  may  give  rise  to 
internal  forces  which  bring  about  failure  of  the  bonds  between  the  con¬ 
stituents,  leading  to  significant  degradation  of  the  strength  and  use¬ 
fulness  of  the  parj:  or  structure. 

In  addition  t?  melting,  vaporization,  and  fracture  due  to  thermal 

I* 

stress,  damage  may  also  be  brought  about  by  other  means.  If  the  heated 
object  is  a  load  carrying  member,  the  interaction  of  the  thermal  effect 
with  the  load  must  be  considered.  The  increase  in  temperature  may  cause 
such  a  reduction  in  tho  yield  strength  of  the  material  that  the  load 
carrying  ability  of  the  part  or  structure  nay  be  seriously  reduced,  causing 
failure  to  occur  at  a  load  which  would  be  in  the  safe  range  for  the  un¬ 
heated  structure.  Creep  rates  may  also  be  higher  at  the  elevated  temp¬ 
eratures  brought  about  by  laser  heating.  Damage  may  also  ba  caused  in 
an  indirect  menner,  as  for  example,  if  the  laser  heating  brings  about  the 
removal  of  a  protective  coating. 
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In  view  of  these  many  possible  consequences  of  heating  of  solids, 
it  is  imperative  that  methods  for  determining  the  resulting  temperature 
distributions  be  available.  For  a  material  which  is  thermally  isotropic 
and  homogeneous,  so  that  conductivity,  density  and  specific  heat  are  each 
a  single  constant,  the  temperature,  as  a  function  of  time  and  position, 
is  the  solution  of  the  familiar  heat  conduction  equation.  If,  in  addition, 
all  thermal  properties  are  independent  of  temperature,  the  equation  is 
linear  with  constant  coefficients.  Solutions  for  some  interesting  and 
important  problems  may  be  obtained  in  closed  form,  or  in  terms  of  a 
Fourier  series.  Many  such  solutions  have  been  collected  by  Carslaw  and 
Jaeger  j[l].  On  the  other  hand,  if  the  thermal  properties  vary  from  point 
to  point,  or  are  functions  of  temperature,  the  possibility  of  finding  such 
solutions  becomes  very  unlikely. 

Problems  for  which  an  exact  .solution  cannot  be  obtained  may  be  attacked 
by  approximate  or  numerical  methods.  Finite  difference  methods  have  been 
used  extensively,  and  have  been  applied  to  both  one  and  two  dimensional 
problems  (.2-6] .  In  addition,  finite  element  methods  for  heat  conduction 
calculations  have  been  developed  in  the  past  decade  [7-9].  Surveys  of 
recent  work  are  also  available  [10,11].  One  disadvantage  of  the  numerical 
method  is  that  thermal  stresses  cannot  be  determined  easily  with  good 
accuracy,  for  the  thermal  stresses  depend  on  the  gradients  of  temperatures, 
and  these  cannot  be  determined  with  good  accuracy  by  a  method  designed 
to  yield  numerical  values  for  the  temperatures  themselves.  A  variational 
method  for  heat  transfer  problems  has  been  given  by  Biot  [12]  and  has 
baen  applied  to  problems  including  phase  changes  [13,14].  A  closely  re¬ 
lated  approximate  method  has  been  given  by  Goodman  [15,16].  To  date, 
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those  nethods  have  been  applied  to  only  one  dimensional  problems,  but 
could  as  well  be  used  in  nulti-dimensional  heat  conduction. 

The  phenomenon  of  phase  change  is  not  described  by  the  heat  con¬ 
duction  equation,  per  se,  but  problems  involving  a  change  in  phase  have 
been  treated  by  introducing  a  moving  boundary  separating  the  two  phases, 
and  attempting  to  solve  the  heat  conduction  equation  in  each  phase,  with 
the  appropriate  matching  being  enforced  at  the  interface.  Problems  in¬ 
volving  such  a  moving  boundary  are  inherently  non-linear,  and  plane  one 
dimensional  problems  have  been  treated  most  extensively  [3,  13,  14,  15, 

16]  by  the  above  mentioned  methods,  and  others  [17-20].  In  this  case,  the 
moving  boundary  is  plane.  In  two  or  three  dimensional  heat  conduction, 
both  the  location  and  the  shape  of  the  boundary  are  unknown,  leading  to 
much  more  difficult  problems  [21,22].  Multiple  phase  changes  have  also 
been  considered  [25]. 

In  real  materials,  the  thermal  properties  are  functions  of  temperature, 
and  it  is  unlikely  that  any  exact  solutions  will  be  found,  although  the 
approximate  methods  can  be  applied  [24],  The  possibility  of  finding 
exact,  or  even  good  approximate  solutions,  for  heat  conduction  problems 
in  composite  or  layered  materials  is  quite  unlikely.  Such  problems  can 
be  wore  successfully  treated  by  one  of  the  numerical  methods. 

In  applying  the  method  of  finite  differences  one  first  writes  down 
the  governing  differential  equation,  and  then  seeks  to  develop  a  set  of 
difference,  or  algebraic  equations  through  the  process  of  discretization. 

In  the  finite  element  method,  as  used  elsewhere,  [9],  a  simple,  approximate 
temperature  distribution  is  assumed  for  each  finite  element,  and  defined 
in  terms  of  nodal  temperatures.  The  resulting  temperature  field  is  sub¬ 
stituted  into  a  global  energy  functional,  the  first  variation  of  which 


leads  to  a  set  of  algebraic  relationships  between  the  nodal  temperatures ,  i 

l 

at  each  tine.  In  a  transient  problem,  this  system  of  equations  must  be 
solved  at  each  time  step,  necessitating  the  inversion  of  a  large  matrix 

i 

at  each  time  increment.  A  finite  element  approach  is  expected  to  be 
particularly  advantageous  for  heat  conduction  in  an  inhomogeneous  material, 
whether  the  inhomogenieities  are  initially  present,  or  arise  from  temp¬ 
erature  dependent  properties,  for  the  thermal  properties  of  each  element 
need  not  be  the  same,  or  constant.  Hie  possibility  of  treating  a  phase 
change  is  also  apparent,  through  allowing  certain  of  the  finite  elements 
to  be  of  one  phase,  while  other  elements  are  of  another  phase.  In  such  an 
approach,  however,  one  would  wish  to  use  extremely  small  elements  so  that 
the  phase  boundary  can  be  located  with  good  accuracy.  In  the  finite 
element  methods  previously  developed,  the  goal  has  been  to  use  as  few 
elements  as  possible,  because  of  computer  storage  limitations  and  the  long 
running  times  required  to  solve  a  large  system  of  equations  at  each  time 
increment.  Consequently,  in  the  present  work,  a  different  finite  element 
approach  was  used,  one  which  does  not  require  the  solution  of  a  large 
system  of  equations  at  each  time  increment.  For  the  same  number  of  elements, 
the  present  method  is  therefore  expected  to  be  less  accurate  than  the  tra¬ 
ditional  finite  element  methods,  but  is  expected  to  be  more  feasible  for 
the  approximate  determination  of  the  liquid-solid  interface  in  two  dimen¬ 
sional  heat  conduction  problems  with  change  of  phase. 

In  the  section  to  follow,  a  general  method  for  solving  heat  conduction 
problems  with  prescribed  boundary  flux  will  be  given.  The  geometry  to  be 
considered  is  of  practical  interest,  namely,  a  plate,  slab,  or  sheet  of , 
uniform  thickness,  subject  to  a  prescribed  flux  on  one  face.  In  the  third 
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section,  some  comparisons  between  results  obtained  with  the  new  method  : 
are  compared  with  results  obtained  by  other  methods,  both  exact  and  approx¬ 
imate,  Finally,  in  Section  IV,  some  results  for  two-dimensional  heat 
flow  with  change  of  phase  are  given. 

■  Although  numerous,  the  works  cited  here  by  no  means  give  a  complete 
survey  of  the  methods  available  for  heat  conduction  problems,  nor  of  re¬ 
sults  obtained  for  problems  involving  a  change  of  phase.  For  a  more  com¬ 
plete  treatment  (through  1964) ,  the  review  by  Muehlbauer  and  Sunderland 
[2SJ  may  be  consulted.  More  recent  work  is  described  by  Boley  [26] . 
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II.  DEVELOPMENT  OF  FINITE  ELEMENT  METHOD 

? 

We  consider  a  cylindrical  disk,  with  a  prescribed  flux  having  a 
known  variation  in  the  radial  direction  and  constant  in  tine  applied  over 
the  central  region.  Exact  solutions  to  this  problem  have  been  given  for 
temperatures  below  melting  in  a  homogeneous  and  thermally  isotropic  disk 
[27].  For  sufficiently  large  flux,  the  temperatures  within  the  disk  will 
be  brought  to  (the  melting  point.  The  determination  of  the  temperature 
and  tho  location  of  the  solid-liquid  interface  then  becomes  quite  diffi¬ 
cult,  and  has  been  extensively  considered  only  for  one  dimensional  prob¬ 
lem,  such  as  occur  when  an  infinite  slab  is  subject  to  a  uniform  flux*.1 

The.  approach  to  be  taken  is  to  divide  the  region  of  interest  into  a 
large  number  of  finite  segments  and  .perform  a  heat  balance  on  each.  We 

‘  i 

begin  by  dividing  the  thickness,  l,  of  the  disk  into  a  number,  N,; ,  of 
layers  which,  for  convenience,  will  all  be  taken  to  be  of'  the  same  thick¬ 
ness,  although  each  layer  can  easily  be  assigned  a  different  thickness,  if 

desired.  In  a  similar  manner,  the  radius,  a,  is  divided  into  Ng  equal  seg- 

« 

stents,  although  unequal  segments  may  again  be  used,  if  desired.  For  the  > 
case  of  an  axially  symmetric  flux,  there  will  be  no  circumferential  flow 
of  heat,  nor  any  circumferential  temperature  variations,  and  consequently, 
no  subdivision  in  the  circumferential  direction  is  required.  Hence,  the 
disk  has  been  divided  into  Ng  x  NR  segments,  each  of  which  isa  in  the  form 
of  a  solid  annulus.  In  order  to  provide  for  eventual  inhomogeneities  in 

the  thermal  properties .  which  may  arise  either  through  consideration  of 

| 

a  composite  material  or  as  a  consequence  of  material  properties  which  are 

i 

a  function  of  temperature,  or  because  of  a  phase  change,  it  is  found  to 
be  convenient  to  define  certain  arrays  (Ngx  Nq)  describing  the  properties 
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of  each  finite  element.  In  ail  that  As  to  follow,  the  first  index  de¬ 
notes  the  row  number,  as  in  Figure  1,  and  is  related  to  the  %  or  thickness 
coordinate,  and  the  second  index  denotes  the  coluan  number,  and  is  related 
to  the  radial  or  r  coordinate.  Hie  array  shown  in  Figure  1  nay  be  thought 
of  as  the  right  half  of  a  diametral  section  of  the  disk. 


r^  •  inner  radius  of  cells  in  coluan  j,  cm 

•  outer  radius  of  cells  in  column  j,  css. 

..  Az  ■  area  of  interface  separating  cells  in  coluan  j 

*  2  2  2  1 

■  ),  ca  . 

°j.  4j 

A-  «  Area  of  interface  separating  cells  in  coluan  j  from  cells  in 
coluan  j  *  1 

•  2«r0^*/NR,  ca2 

p.j  •  density  of  cell  of  row  nuaber  i  and  coluan  j,  gm/cm3 


M.j  *  mass  of  cell  i,  j,  ga. 


*  p  Vh 
iizi 

*  heat  of  fusion  for  cell  i,  j,  joules/gm 

K^,  *  conductivity  of  cell  i,  j  in  radial  and  axial  directions 

respectively,  joules/ (sec  cm  *C) 

Qn  ,  Q-  *  rates  of  heat  flow  out  of  cell  i,  j  in  radial  and  axial 

TJ  hj 

directions,  respectively,  joules/sec. 

©ij  *  Average  temperature  in  cell  i,  j,  *C  ' 

Cp^  >  specific  heat  of  cell  i,  j,  joule/ (ga*C). 
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Thu  heat  fluxes  for  a  typical  element  are  shown  graphically  in  Figure  2. 

* 

The  rate  at  which  heat  flows  radially  from  cell  i,  j  into  cell  i,  j+1  is 
given  approximately  by  the  product  of  the  conductivity,  the  area  of  the 
interface,  and  an  approximate  value  for  the  temperature  gradient  at  the 
interface  between  these  two  cells,  or 


Qr  «  Kr.  Ar  I  jjd-L  1 

ij  Rij  Rj  l  a/Nc  J 


(1) 


Similarly,  the  rate  at  which  heat  flows  from  the  cell  i,  j  axially  into 
the  cell  i+1,  j  is  approximately 


Q?  * 

1  *  *1 


l  t/NR  5 


(2) 


‘ij  *ij  “i 

In  addition  to  the  flow  of  heat  from  one  cell  to  another  through  conduc¬ 
tion,  the  heat  balance  must  take  into  account  the  possibility  of  heat 
addition  to  each  cell  from  external  sources,  and  from  internal  reactions, 
through  a  heat  addition  term  for  each  cell  which  includes  both  effects. 


Fij  *  FE 


♦  Ft 


(3) 


ij  *3 

Fjj ,  the  heat  added,  has  dimensions  of  joules/sec.  The  heat  balance  is 
performed  for  each  cell  by  adding  the  heat  entering  cell  i,  j  from  the  neigh¬ 
boring  cells  i-1,  j  and  i,  j-1  to  the  flux  added,  Fjj,  and  subtracting 
the  flux  leaving,  Qr^  and  Qz^ .  If  cell  i,  j  is  at.  a  temperature  lers 
than  the  phase  transformation  temperature,  the  temperature  will  rise  in 
a  time  At  an  amount  given  by 

A^ij  *  AQ*j  *  ^t/(Cp^  (4) 

■  Fij  *  s.j..  *  -  V, .  (5i 

The  temperature  rise  will  continue  in  later  time  increments  At  un¬ 
til  such  time  as  tho  temperature  at  which  the  phase  transformation  takes 
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place  is  reached.  The  total  energy  which  must  be  added  after  the  melting 
temperature  is  reached  in  order  to  bring  about  a  complete  change  of  phase 
of  all  the  mass  in  cell  i,  j  is  given  by  It  is  convenient  to 

compute  this  quantity  for  each  cell,  and  store  the  result  as  an  array  Py . 
Then,  beginning  with  the  first  time  increment  in  which  the  temperature  of 
the  cell  i,  j  is  at  or  above  the  melting  temperature  ^m»  the  stored  quan¬ 
tity  Pj^  is  reduced  from  its  current  value  by  the  amount  AQ^At  in  the 
time  increment  At  and  the  temperature  is  held  fixed  at  the  melting  temp¬ 
erature.  At  K  time  increments  after  cell  i,  j  first  reaches  T#, 

pij  "  pii  -  k2  j  ">ljk  “k  <« 

and  the  entire  cell  is  presumed  to  complete  the  phase  change  during  that 
time  increment  in  which  *►  o.  The  material  in  cell  i,  j  is  then  per¬ 
mitted  to  increase  in  temperature  once  again,  with  specific  heat  and  con¬ 
ductivity  now  being  appropriate  for  the  new  phase,  if  it  is  assumed  to 
remain  on  the  surface.  Alternatively,  the  material  which  has  undergone 
the  phase  change  may  be  assumed  to  be  instantaneously  removed,  as  in  the 
case  of  ablation,  or  melting  where  sufficient  forces  are  present  to  over¬ 
come  the  surface  tension  of  the  melt  and  remove  the  melted  material  as 


it  is  formed.  In  the  case  of  complete  melt  removal  (or  equivalently, 
ablation),  as. As  considered  in  all  of  the  examples  to  be  discussed  in  the 
next  sections,  the  conductivities,  and  Kz^  are  set  to  zero  at  the 
time  the  material  is  removed,  which  has  the  effect  of  removing  the  cell 


t,  j  from  the  heat  balance  calculations. 


The  flux 


» 


representing  the 


flux  added  externally  to  cell  i,  j  is  transferred  to  the  cell  immediately 


below,  i.e, ,  to  cell  i-tl,  j  at  the  same  time. 


n 


,.y>  'tnitn,' 


The  heat  balance  for  cells  or  elenents  on  the  boundary  must  be  handled 
somewhat  differently  from  those  on  the  interior  so  as  to  provide  for  the 
approximate  satisfaction  of  any  given  boundary  conditions.  In  the  case 
of  a  prescribed  boundary  temperature,  the  temperature  of  that  boundary  cell 

•  i 

is  held  fixed  for  all  time  increments.  In  the  case  of  an  insulated  boundary, 
the  net  flux  in  the  direction  normal  to  the  boundary  is  set  to  zero.  In 
the  case  of  a  prescribed  flux,  the  net  flux  in  the  direction  normal  to 
the  boundary  is  set  to  the  desired  value,  and  in  the  case  of  a  boundary 
condition  described  by  Newton's  law  of  cooling,  the  flux  normal  to  the 
boundary  is  set  to  the  difference  between  the  temperature  of  the  boundary 
cell  and  the  surroundings,  with  the  appropriate  multiplicative  constant 
included.  By  symmetry,  the  radial  flux  entering  cells  in  column  1  is 
zero. 

Because  of  the  approximate  way  in  which  the  conductibn  is  computed 
(Equations  1  and  2)  the  time  increment  At  must  be  sufficiently  small  in 
order  to  insure  stability  of  the  solution.  The  flux  for  the  entire  time 
increment  is  computed  based  on  the  temperature  difference  between  adjacent 
cells  at  the  beginning  of  the  increment.  It  is  evident  that  if  the  time 
increment  is  permitted  to  be  too  large,  sufficient  heat  would  flow  to  cause 
the  initially  colder  element  to  become  hotter  than  the  initially  hotter 
element,  which  is  obviously  impossible.  A  physical  derivation  of  the 
appropriate  stability  condition  is  easily  obtained.  If  two  cells,  the 
centers  of  which  are  a  distance  Az  apart  are  initially  at  temperatures 
Tj  and  Y2,  and  if  the  conductivity  is  K,  in  time  At  the  total  heat  flow 
will  be 


Q»  KA  At  (Tj  -  T2)/Az 


(7) 


where  A  is  the  area  pf  the  annulus  separating  the  two  eleaents.  If  the 

specific  heat  (per  unit  volume)  of  each  is  pC  ,  then  at  the  above  flux, 

P 

temperature  equilibrium  between  the  two  elements  is  reached  when 

Q  ■  pCpV  (Tj  -  T2)/2  (8) 

where  V  is  the  voluae  of  each  element,  approximately  AAz  for  each*  Elim¬ 
inating  Q  between  Equations  7  and  8,  we  find 

At  -  pCp  (As)2/(2K)  (9) 

But  temperature  equilibrium  cannot  be  established  in  a  finite  time,  hence 
At  must  be  less  than  this  bound.  Moreover,  in  two  dimensional  heat  flow, 
both  components  of  flux  mist  be  considered,  leading  to  a  more  stringent 
requirement.  However,  if  the  two  element  spacing  distances,  a/Nc  and 
i/Ng  are  not  comparable,  the  maxima  permissible  time  increment  is  sat¬ 
isfactorily  determined  by  using  the  lesser  distance  increment  in  Equation  9. 

* 

A  method  given  by  Dusinberre  (28]  for  the  hand  calculation  of  transient 
temperatures  is  very  similar  to  the  method  used  here  to  compute  the  heat 
conduction.  If  the  thermal  properties  are  uniform  and  constant,  the  method 
can  also  be  shown  to  generate  the  same  equations  as  would  arise  with  the 
method  of  finite  differences,  with  central  differences  being  used  in  the 
space  variables  and  a  forward  difference  employed  in  time.  The  stability 
und  convergence  of  such  techniques  applied  to  linear  problems  (constant 
properties  and  no  melting)  have  been  extensively  studied  [29].  Numerical 
methods  for  the  treatment  of  the  non-linear  problems  have  also  been  dis¬ 
cussed  (30J . 

A  computer  program  was  written  and  used  to  perform  the  calculations 
described  in  the  next  section.  A  listing  of  that  program,  together  with 
a  brief  documentation,  is  given  as  an  appendix. 
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III.  COMPARISONS  WITH  KNOWN  SOLUTIONS 


We  will  now  consider  several  example  problems,  making  comparisons 
between  temperature  fields  predicted  by  the  finite  element  method  described 
in  the  preceeding  section,  and  exact  solutions,  when  possible,  and  other 
approximate  methods,  when  available. 

a.  One  Dimensional  Heating  of  Thick  Slab 

As  a  first  example,  wo  will  consider  one  dimensional  heating  of  a 
thick  (infinite)  slab  having  properties  approximating  those  of  A1203. 

We  take 

p  «  3.8  gm/ca3 

K  ■  0.104  joule/ (cm  sec  *C)  ! 

Cp  «  0.88S  joule/ (gm*C) 

Tn  «  2313*K 

and  consider  a  uniform  flux  of  4000  joules/ (cm2sec)  striking  material 
initially  at  300*K.  The  exact  solution,  valid  to  the  onset  of  melting, 
is  known  [1]  to  be 

1/2 

T.T0.SS2_i„fc{_ij7r}  CW) 

where  a  is  the.  diffusivity,  K/ (pCp) ,  and  ierfc  denotes  the  first  integral 
of  the  complementary  error  function.  The  front  surface  temperature  is 
given  by 

T  -  T0  -  IE  (at/t)1/2  (11) 

end,  for  the  parameters  given  above,  the  surface  melts  at  0.0695  sec. 

The  temperature  profiles  given  by  Equation  10  are  plotted  at  0.02, 
0.04,  and  0.06  seconds  in  Figure  3  as  the  solid  lines.  Temperature  pro¬ 
files  for  the  same  heating  conditions  were  determined  by  the  finite  ^ele¬ 
ment  method  through  considering  a  slab  0.2  cm  thick  (from  Figure  3,  it 


can  be  seen  that  a  slab  of  such  thickness  nay  be  regarded  as  being  of 

infinite  thickness  for  tines  less  than  0.06  seconds).  The  slab  was 

! 

divided  into  20  and  40  layers,. with  tine  steps  of  0.001  and  0.00025  seconds 
used  in  the  two  cases,  respectively.  (Tine  steps  less  than  0.0016  and 
0.0004  are  required  here  for  stability).  The  average  tenperature  in  the 
.  first  0.01  cn  of  the  naterial  at  t  ■  0.06  was  found  to  differ  by  less  than 
f  0.2V  between  the  two  calculations.  The  results  of  the  finite  elenent  ’< 

calculations  using  the  larger  layer  thickness  and  large  tine  step  are 

* 

given  as  the  circles  in  Figure  3.  As  can  be  seen,  the  results  of  the 
.  exact  and  approximate  calculations  are,  for  the  nost  part,  indistinguishable, 
b.  One  Dimensional  Melting  of  Finite  Slab 

As  a  second  example,  we  consider  the  one-dimensional  heating  and  change 

of  phase  of  a  finite  slab,  assuning  the  nelt  is  instantaneously  removed* 

For  a  prescribed  uniform  flux  F  on  the  surface  x  «  i,  and'  the  other  face 

(x  •  0)  insulated,  the  resulting  one  dimensional  temperature  field  is  [1] 

T-T  *  lUl  *  f  .  Jj,  i  e  “an^w  cos  (12a) 

0  pCpt  K  l  TP  w*  n»l  n2  '  1  } 

.  EdbE  ierfe  [£2i!2£2]  l  (12b) 

K  «■»  2/oF  2 /St1 

valid  until  such  time  as  melting  begins  at  the  front  surface.  Masters  [3] 

has  given  results  obtained  by  a  finite  difference  calculation  for  the 

predicted  recession  of  the  front  surface  for  an  aluninun  slab,  using  the 

properties 

Tm  »  993*K 

L  •  418  joules/gn 

X  •  2.09  joul*/(cn  sec  *C) 
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a  •  K/(pCp)  *  1  ca2/sec 
t  *  0.3  cm 

We  assume  the  complete  and  instantaneous  renoval  of  the  Melt,  resulting 
from  the  application  of  a  flux  of  41,800  joules/ (cm  sec).  Taking  the 
density  to  be  2.7  ga/ca2,  the  value  of  specific  heat  required  to  yield 
«  unit  diffusivity  is  C..  *  0.773  Joule/  (gn*C) 

The  total  tine  required  to  aelt  a  slab  under  one  dimensional  heating 
is  readily  determined  from  an  overall  heat  balance  to  be 

Ft  -  (L  ♦  Cp(TB  -  T0)Jpi  (13) 

or  0.0186  seconds.  From  the  exact  solution,  (Equation  12),  the  tine  re¬ 
quired  for  the  front  surface  to  reach  the  melting  solution  was  found  to 
be  0.0011  sec.,  hence  0.0175  sec.  are  required  to  move  the  melting  line 
froa  the  front  to  the  rear  surface.  This  tiae  was  also  computed  by  the 
finite  eloaent  method,  first  by  dividing  the  slab  into  10  layers  and  using 
a  tiae  increment  of  .200  u  sec  and  by  dividing  the  slab  into  20  layers  and 
using  a  tiae  increment  of  SO  u  sec.  The  two  aelting  tines  were  determined 
to  be  0.0182  see  and  0.0178  sec,  indicating  that  the  approximate  method 
appears  to  be  converging  to  the  correct  solution  as  the  size  of  the  incre- 

i 

nents  is  reduced.  In  both  of  these  cases,  the  time  steps  used  were  just 
under  50%  of  the  minimum  tiae  step  required  for  stability.  The  location 
of  the  moving  free  surface  is  shown  in  Figure  4.  The  solid  line  is  the 
result  given  by  Masters  [3j,  and  implies  a  predicted  tine  (>  0.02  sec)  which 
is  greater  than  the  tiae  deduced  from  the  overall  heat  balance  (0.0186  sec). 
The  other  two  curves  in  Figure  4  depict  the  results  of  the  finite  eleaent 
calculation,  using  20  layers  and  At  *  50  p  sec.  The  lower  curve  depicts 
the  tiae  at  which  various  points  first  reach  the  aelting  temperature,  and 
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the  upper  curve  indicates  the  tines  at  which  particles  are  completely 
neited  and  removed.  The  interval  between  the  curves  indicates  the  region 
in  the  process  of  melting.  . 

The  problem  of  one  dimensional  melting  with  complete  removal  of  .the 
melted  material  has  also  been  treated  by  Citron,  who  has  given  [14]  an 
approximate  (rather  than  numerical)  method  and  some  typical  results.  His 
results  are  presented  in  such  a  manner  as  to  indicate- that  only  a  single 


dimensionless  parameter,  M/r,  is  required  to  characterize  the  problem. 

M  ..  CnlF  r  «  kT  / 

r  2  XL  *  tP  !  * 


The  value  of  M/r  for  the  parameters  used  in  the  previously  described  cal¬ 
culation  is  9.66.  These  calculations  were  then  repeated  with  a  flux 

t 

F  ■  20,750  joule/ (cm2sec)  so  as  to  give  a  value  of  M/r  ■  4.89  identical  to 
the  one  used  by  Citron.  The  two  dashed  curves  of  Figure  5  are  results  given 
by  Citron  as  first  and  second  approximations.  The  abcissa  is  a  dimension¬ 


less  time,  defined  by 


t  «  (t  -  t*) 

J2 


t*  being  the  time  of  front  surface  melting. 

The  ordinate  is  the  fraction  of  the  thickness  which  has  melted  at 


dimensionless  time  r.  The  solid  curve,  labeled  r'»  0.235,  depicts  the  re¬ 
sults  obtained  by  the  finite  element  method,  and  should  be  in  agreement 
with  the  results  of  Citron.  The  apparent  discrepancy  at  the  early  time 
was  found  to  bo  due  to  the  inapplicability  of  the  Citron  method  at  small 
values  of  r,  even  though  r  does  not  appear  explicitly  in  the  results.  In 
developing  this  approximate  method.  Citron  assumed  the  exponential  terms 
in  Equation  12a  to  be  negligible  at  the  melting  time,  t*.  With  this 
assumption,  we. have 

r  ■ 
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A  r  of  0.235  leads  to  a  meaningless  negative  time  of  first  melting.  Hie 
approximation  made  by  Citron  is  valid  only  for  r  -  2/3  or  greater.  In 
order  to  compare  the  results  of  the  present  method  with  those  given  by 
Citron  with  a  value  of  r  such  that  his  results  are  valid,  calculations 
were  performed  for  a  flux  of  F  *  20,750  joule/ (cm2sec)  and  all  properties 
those  of  .3  cm  thick  aluminum,  except  that  a  melting  temperature  of  2373*K 
was  used.  This  leads  to  a  value  of  r  of  0.705,  with  M/r  •  4.89.  The 
results  of  this  calculation  are  given  in  the  Figure  as  the  indicated 
solid  line.  Although  there  is  some  difference  between  the  two  pre¬ 
dictions,  they  are  now  in  qualitative  agreement. 

•  c.  Two  Dimensional  Heating  of  a  Thin  Sheet 


As  a  third  example,  we  consider  the  two  dimensional  temperature  dis¬ 
tribution  resulting  from  a  steady  axi-symmetric  flux  of  a  Gaussian  distri¬ 


bution  acting  on  a  thin  sheet. 

F(r)  -  Fe 


-r2/2o2 


The  pre-melting  solution  can  be  compared  with  the  temperature  profiles 
predicted  by  a  computer  program  written  at  AFWL  [31]  which  evaluates  the 
solution  given  by  Oicer  [27]. 

As  an  example,  we  will  consider  a  titanium  sheet,  0.04  cm.  thick  and 
5  cm,  in  diameter,  initially  at  300*K,  having  thermal  properties 
p  *  4.43  gm/cm3 
K  ■  0.145  joules/ (cm*Csec) 

Cp  •  0.77  joules/(gm*C) 

Tm  ■  1900*K 
L  ■  390  joules/gm 

We  assume  the  dis\  is  subjected  to  a  beam  of  peak  intern*  \iy  of  2000 
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joules/ (sec  cm2)  with  a  *  0.25  cm.  The  pre-melting  temperature  field  was 
computed  using  the  finite  element  method  by  (a)  dividing  the  slab  into 
10  layers  and  the  radius  into  40  annular  regions,  using  a  time  step  of 
80  u  sec,  and  (b)  by  dividing  the  sheet  into  20  layers  and  20  annular 
rings,  using  a  time  step  of  20ysec.  In  both  cases,  the  time  step  used 
was  about  40%  of  the  minimum  value  necessary  to  satisfy  the  one  dimensional 
stability  criterion.  In  the  most  extreme  case,  AR/Az>15,  and  it  was  assumed 
that  instabilities  would  primarily  arise  due  to  the  axial,  rather  than  the 
radial,  flux.  The  temperature  profiles  at  z  *  0.002  cm  (A/20)  were  found 
to  differ  by  less  than  1%  for  all  r,  at  tines  t  •  0.03  and  t  *  0.10  sec. 
Temperature  profiles  as  computed  with  the  finite  element  method  (solid  lines) 
are  compared  with  the  results  of  the  Fourier  series  solution  (AFNL  computer 

\ 

program)  in  Figure  6.  The  agreement  is  excellent,  particularly  within  the 
beam  radius  (2o  «  0.5  cm).  Temperature  profiles  through  the  slab  thick- 
ness,  at  r  *  0.0625  cm,  are  compared  at  two  values  of  time  in  Table  I. 

The  exact  solution  indicated  that  melting  begins  at  the  front  surface  at 
t  *  0.1037  sec.  The  two  finite  element  solutions  gave  0.1042  and  0.1064 
sec. ,  respectively.  Although  this  is  a  two  dimensional  problem,  axial 
heat  flow  predominates  and  the  one  dimensional  approximation  is  quite 
satisfactory  for  it  leads  to  a  predicted  time  to  first  melting  which  is 
only  10%  low. 
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FINITE  ELEMENT  SOLUTION 


'  IV.  APPLICATIONS 

In  order  to  gain  further  experience  in  the  utilization  of  the  new 
Method,  several  two  diaensional  heat  conduction  and  melting  problems  were 
studied.  Applications  involving  several  materials  of  interest  were  con¬ 
sidered,  viz.  AI2O3,  Titanium,  Stainless  Steel,  and  Magnesium*  In  none 
of  these  cases  is  another  solution  available  for  comparison,  although  ( 
experimental  results  corresponding  to  several  of  these  cases  are  available, 
a.  Melting  of  Thick  Alumina  Slab 

The  predicted  time  of  melting  through  a  thick  (0.952  cm)  slab  of 

AI2O3,  5  cm  in  diameter,  due  to  an  absorbed  flux  of  peak  intensity  or 

2 

4000  joules/ (cm  sec)  and  Gaussian  parameter  of  0  ■  0.70  cm  was  computed 

with  a  number  of  runs.  In  all  cases,  the  following  parameters  were  assumed 
3 

p  «  3.8  gm/cm 

L  >  1070  joules/gm 

K  *  0.104  joules/ (cm  sec  #C) 

Cp  •  0.885  joules/ (gm  *C) 

T#  «  2313  *K 

and  the  initial  temperature  was  taken  to  be  300*K.  The  slab  was  divided 
into  various  (10,20,  and  40)  layers,  and  the  radius  into  20  and  40  segments. 
The  time  steps  which  satisfy  the  one  dimensional  stability  criterion  are 
.147,  .0368  and  .0092  sec.,  respectively.  The  results  of  the  various  compu¬ 
tations  are  given  as  Table  II.  The  temperature  at  r»a  and  z«0  at  the  time 
of  complete  melt  through  was  found  to  have  increased  only  2*C,  hence  these 

I 

results  are  also  applicable  to  disks  of  any  larger  diameter.  Excellent 

1 

agreement  between  the  results  of  various  runs  is  evident,  despite  the  sig- 
nifeiant  variations  in  the  size  of  space  and  time  increments.  Instantaneous 
removal  of  the  melted  material  was  assumed. 
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Table  II 


Tino  Required  to  Melt  A^Oj  Slab,  0.953  cm  Thick  end  2*5  cm  DiaMeter 

2 

with  F  *  4Kw/cn  and  a  *  0.70  cm 


Run  )io. 

Nunbor  Layers 

Number  Radial 
Segments 

Tiee  Step 
Sec. 

Melt  Through 
Time  Sec. 

1 

10 

20 

*  .01 

2.86 

2 

10 

40 

.01 

2.83 

3 

20 

20 

.005 

2.845 

4 

20 

40 

.005 

2.85 

5 

40 

40 

.0025 

2.8475 
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b.  Melting  of  a  Thin  Titaniun  Sheet 

We  consider  now  the  consequence  of  a  Gaussian  bean  of  t  •  0.25  being 
absorbed  by  a  0.04  cm  thick  sheet  of  Titaniun,  initially  at  300*K.  We 
tako, the  thormal  properties  to  be  as  in  Section  IIlc  and  At  ■  80  u  sec. 

Tho  times  at  which  the  various  finite  segments  were  found  to  be  completely 
nelted  (and  assumed  to  be  removed)  are  depicted  graphically  in  Figure  7. 

The  tines  given  in  the  Figure  are  measured  fron  the  onset  of  melting  at 
the  front  surface.  The  exact  and  discretized  approximation  to  the  prescribed 
flux  is  also  shown.  The  propagation  of  the  free  boundary  through  the  sheet 
is  shown  as  the  solid  lines  in  Figure  8  for  peak  absorbed  intensities  of 
1500,  2000,  and  2500  joules/ (cm2sec).  For  purposes  of  comparison,  the  pre¬ 
dicted  melting  rate  for  one  dimensional  heating  with  a  peak  intensity  of 
2000  joules/ (cm2soc)  is  also  given  as  the  dashed  line.  As  was  noted  in 
Section  IIlc  a  Titanium  sheet  of  this  thickness  can  be  quite  well 
approximated  by  a  one  dimensional  problem,  despite  the  two  dimensionality 
of  the  temperature  distribution,  evident  from  Figure  6.  Although  no 
analytical  results  for  this  problem  are  available,  the  predicted  melt 
through  time  is  in  satisfactory  agreement  with  experimental  results  [31]. 

c.  Melting  of  Thin  Stainless  Steel 

In  an  experimental  test  program  conducted  at  WPaFB,  a  16  mil  sheet  of 
304  Staialess  Steel  was  irradiated  by  a  beam  having  a  total  power  of  9 
kilowatts  and  a  diameter  of  2.44  cm.  Temperatures  wer6  measured  by  a 
thermocouple  attached  to  the  rear  surface.  The  results  [31]  were  as 

j 

indicated  by  the  circles  in  Figure  9.  From  notion  picture  films  taken  of 
the  event,  the  tine  at  which  complete  melt  through  occured  was  estimated 
to  be  0.40  sec.  Assuming  a  Gaussian  bean  profile,  a  peak  intensity  of 
3850  joules/ (cm2sec)  was  estimated  from  the  relationship 


27 


018  .053  .172 


038  8S0*  JLV  0o9~lV 


SSo2cMro^io<o 

CM  cvl  OJ  CM  CVS  CVI  CM 


fO  o 
6  * 


~»?0100(0010 
<D<DNCDC0(7>0>qO 
OOOOOOO  —  — 


v  i/> 

06  M 


VI 

i  : 


o  u>  o>  £  —  ioa> 

{\}rororo^Ju>u>»o 

ooooooooo 
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The  material  was  assumed  to  be  removed  immediately  after  melting. 

Calculations  of  the  temperature  distribution  and  progression  of  melting 

were  made,  using  the  following  thermal  properties: 

p  *  7.9  ga/crn3 

K  *  0.24  joules/(sec  cm  #C) 

C  •  0.42  joules/ (gm*C) 

P 

L  *  390  joules/gm 
Tm  •  1700#K 

An  initial  temperature  of  300*K  was  assumed  and  a  disk  radius  of  5  cm  was 
used.  The  disk  was  divided  into  10  layers  and  20  annular  rings,  and  a 
time  step  of  SOusec  was  used  in  the  calculations.  Absorbtances  of  0.2  and 
0.3  were  first  assumed,  leading  to  the  predicted  rear  surface  temperatures 
shown  as  solid  lines  in  Figure  9.  The  predicted  melt  through  times  for 
these  absorbtances  were  0.387  sec  and  0.2S6  sec,  respectively.  While  an 
absorbtonce  of  0.2  leads  to  a  melting  time  which  is  close  to  the  observed 
value,  the  predicted  temperature  history  at  the  rear  surface  is  significantly 
different  from  that  which  was  observed.  The  observations  appear  to  correspond 
to  an  absorbtance  of  about  0,0.8  at  early  times,  followed  by  a  transition 
to  a  higher  value,  perhaps  0,24.  Based  on  this  observation,  the  computer 
program  was  modified  so  that  the  flux  entering  a  cell  on  the  surface  would 
increase  from  the  lower  to  the  higher  value  when  the  average  temperature 
of  that  cell  reached  6Q0*K.  The  resulting  predicted  time  for  complete 
melt  through  is  0,416  sec.,  and  the  thermal  history  of  the. rear  surface  can 
be  seen  from  the  dashed  curve  to  be  in  satisfactory  agreement  with  the 
experimentally  determined  results.  These  results  would  seem  to  indicate 
that  a  significant  change  in  absorbtance  of  Stainless  Steel  takes  piece 
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at  temperatures  well  below  the  melting  temperature.  It  is  conjectured 
that  this  change  takes  place  due  to  the  formation  of  an  oxide  layer  on  the 
heated  surface.  The  formation  of  such  an  oxide  layer  can  be  expected  to 
depend  not  only  on  the  current  temperature,  but  also  on  the  temperature 
history  of  the  surface,  which  would  be  dependent  on  the  magnitude  of  the 
applied  flux. 

d.  Melting  of  Magnesium  Sheets 

Since  the  time  required  for  complete  melt  through  of  a  sheet  is  very 
easy  to  calculate  (Equation  13)  if  the  heat  flow  is  entirely  axial,  and 
very  difficult  if  it  is  not,  it  would  be  very  advantageous  if  there  were 
some  means  of  knowing,  a  priori,  whether  or  not  a  one  dimensional 
approximation  is  appropriate  for  any  given  heating  situation. 

A  computational  study,  considering  magnesium  of  thicknesses  from 
0.08  cm  to  1.28  cm,  beam  parameters  (o)  from  .32  and  1.28  cm  and  peak 
intensities  of  1,  3,  S,  and  lOKjoules/Ccm^ec)  was  undertaken.  A  cylindri¬ 
cal  disk,  10  cm.  in  diameter,  and  having  the  following  thermal  properties 
was  considered,  and  the  melt  was  assumed  to  be  instantaneously  removed 

j  1 

p  »  1.77  gm/cm 

| 

K  •  0.96  joules/ (cm  sec  *C) 

Cp  *  1.04  joule/ (gm  °C) 

L  «  338  joules/gm 
«  905*K 

Each  disk  was  subdivided  into  10  layers,  regardless  of  thickness,  and  the 
radius  was  divided  into  20  segments,  except  in  the  case  of  e  *  0.32i  The 
discretisation  of  the  incident  flux  was  found  to  lead  to  significant;  errors 
when  the  parameter  o  is  comparable  to  the  size  of  the  radial  segment',  or 
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greater,  necessitating  the  use  of  40  segments  for  the  case  of  o  *  0.32  cm. 

In  all  cases,  tiae  steps  wore  kept  below  SOI  of  the  ainiaua  value  required 
to  insure  stability. 

The  results  of  these  calculations  are  given  as  Figure  10,  where  the 
tiae  to  aelt  is  given  as  a  function,  of  thickness  for  various  beaa  paraaeters. 
As  can  be  seen  froa  the  figure*  large  spot  site  and  higher  levels  of  flux 
both  lead  to  aelting  tiaes  in  closer  agreeaent  with  the  aelting  tiae  for 
one  di&ensional  problems,  as  given  by  Equation  13. 

it  might  be  expected  that  the  ratio  o/Jt,  being  the  ratio  of  the  only 
two  lengths  present,  would  characterize  the  "degree  of  one-dimensionality" 
of  a  given  heating  situation  for  a  given  material.  This  was  not  found 

to  be  the  case.  The  ratio,  for  example,  of  t/tj  for  F  ■  S  Kw/cm  can  be 

.  •  ’  » 

seen  to  vary  from  1.348  at  o  •  t  ■  0.32  down  to  less  than  1.1  for  o  ■  Jt  • 

» 

1.28.  A  more  promising  aeans  of  interpretation  appears  to' be  a  plot  such 

*  • 

as  is  given  in  Figure  11,  wherein  the  normalized  melting  time  is  given  as 

i 

a  function  of  total  power,  for  various  sheet  thicknesses.  The  tine  to 
melt  through  has  been  normalized  by  division  by  the  time  required  for  a 
uniform  beam  of  the  same  peak  intensity  to  melt  through  a  similar  slab, 
assuming  one  dimensional  heat  flow.  Computed  points  are  shown,  for  the  various 
thicknesses,  as  circles,  squares  and  triangles,  but  all  the  computed 
points  do  not  appear  to  fall  on  smooth  curves.  For  this  reason,  deduc¬ 
tions  about  the  one-dimensionality  of  problems  well  outside  the  range  of 
parameters  considered  in  this  study  should  not  be  attempted,  but  these 
results  do  appear  to  suggest  that  the  total  powex^  is  a  more  significant 
parameter  than  the  peak  intensity,,  spot  size,  or  the  ratio  of  spot  size  1 
to  thickness*  Froa  Figure  11,  it  can  be  seen  that  the  curves  of  t/tj,  , 
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t-TIME  FO 


for  Various  Bom  Parsooters 


vs.  power  for  different  thicknesses  are  all  parallel ,  that  is,  if  each  is 
translated  along  the  horizontal  axis,  the  three  curves  can  be  made  to 
concidc.  Moreover,  the  distance  that  each  curve  Bust  be  translated  is 
inversely  proportional  to  t,  i.e.,  t/tj  ■  f(P/l).  Thus,  rather  than  a 
dependence  of  Belting  time  on  the  various  psraaeters  having  the  fora 

V"  £(°*M)  (i8) 

we  can  write 

tn  -  tjfCP/i)  (19)' 

where  tm  is  the  tine  required  to  aelt  a  seal 1  hole  with  a  bean  of 
intensity  F,  gaussian  parameter  a  in  a  thickness  t  of  some 
material 

tj  is  the  time  required  for  one  dimensional  Belting  of  the  saae 
material,  with  flux  P  and  thickness  t, 

tx  -  tL  *  C_CT.  -  T„)]  at  '  (20) 

•  r 

P  is  the  total  power 

1  P  -  2HF02  (21) 

i  is  the  thickness 

and  f  is  a  single  valued  function,  which  we  may  expect  to  be  differ¬ 

ent  for  each  material. 


iu i'  magnesium,  f (40, 000)  -  1.1  so  the  time  required  for  aelting  can  be  , 

approximated  to  within  10%  by  the  tine  required  for  one  dimension  of  melting  wheneve 
>  40,000  joules/ (cm  sec).  For  comparison,  f(4S00)  -  2.0  indicating 
that  melting  will  required  twice  as  long  as  in  the  one  dimensional  care. 

In  this  case,  50%  of  the  flux  applied  at  the  center  of  the  heated  area 
is  conducted  axially,  and  50%  radially. 
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Ilie  results  presented  in  Figures  10  and  11  should  be  regarded  as 
being  somewhat  tentative,  as  multiple  computations  of  each  case,  using 
various  spatial  and  temporal  step  sizes,  were  not  undertaken  to  insure 
convergence  of  these  numerical  results. 

l 
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V.  CONCLUSION 

In  the  preceding  sections,  a  method  for  determining  numerically  the 
temperature  distribution  within  a  disk  of  arbitrary  dimensions  and  thermal 
properties  subject  to  a  prescribed  flux  over  one  face  has  been  given. 

The  method  takes  into  account  a  change  of  phase,  and  can  be  used  to  predict 
the  time  required  for  partial  or  complete  melting  of  the  material.  The 
method  was  applied  to  a  number  of  examples  and  was  found  to  give  satis¬ 
factory  agreement  with  other  solutions  and  with  experimental  evidence 
when  available. 

Although  no  calculations  for  composite  materials  were  included  in 
the  examples,  tho  method  was  developed  with  such  calculations  in  mind, 
specifically,  so  that  one  or  more  layers  could  be  taken  to  have  the  thermal 
properties  of  a  material  which  might  be  used  as  a  protective  coating. 

The  method  can  also  be  expected  to  handle  a  number  of  other  problems, 
and  can  be  modified  to  incorporate  other  factors.  A  single  phase  change 
from  solid  to  vapor  (rather  than  liquid)  can  be  handled  by  the  method 
without  any  modification.  Two  phase  changes,  as  first  a  transition  from 
solid  to  liquid,  followed  by  heating  of  the  retained  liquid  and  the 
eventual  vaporization  of  the  liquid  can  be  readily  incorporated  into  the 
»s<n"  od,  Provision  was  also  made  for  the  addition  of  heat  to  each  cell,  so 
that  lunt  releasing  reactions  within  the  material  might  be  considered. 

A  time  dependent  prescribed  flux  can  be  handled  without  difficulty,  as 
van  boundary  conditions  other  than  the  completely  insulated  boundary 
which  was  assumed  in  the  examples  considered  here.  Changes  in  flux  as 
a  function  of  surface  temperature,  such  as  can  occur  through  a  temperature 
dependent  absorb tance  can  be  treated  by  this  method,  as  was  demonstrated 
4.0  one  example.  Prescribed  fluxes  which  are  not  axially  symmetric  could 
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in  principle  bo  treated  by  a  finite  element  method  such  as  this,  but!  it  is 
anticipated  that  the  large  number  of  elements  which  would  be  necessary  to 
handle  a  third  component  in  the  heat  flux  vector  would  require  inordinately 
large  computer  storage  and  running  times  for  these  transient  applications. 
Certain  quasi-three  dimensional  problems,  in  which  the  prescribed  flux  i3 
a  two  dimensional  function  of  the  surface  coordinates,  but  for  which  the 
heal  flow  is  primarily  one  dimensional  are  probably  amenable,  however* 

Such  problems  might  ariso  from  the  spatial  and  temporal  fluctuation  of 
a  distribution  of  radiation. 

Significant  remaining  problems  in  determining  the  consequences 
of  laser  heating  of  materials  would  appear  to  be  in  the  development  of 
models  for  predicting  the  rate  of  mass  removal  and  for  predicting  the 
rate  at  which  heat  As  liberaf  d  in  such  materials  as  titanium  and  magnesium, 
as  has  been  observed,  and  in  determining  the  absorbtance  as  a  function  of 
surface  conditions. 
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APPENDIX 


A  Fortran  Program  was  written  to  perform  the  two  dimensional  heating 
and  melting  calculations  by  the  method  described  in  Section  II.  A  brief 
documentation  follows: 


Card  Numbers 
1S-20 


25-45 

48-56 


60 

66-72 

76-82 

85-94 

95-106 


Computation  Perforaed 


Input  data  read  in.  All  temperatures  Kelvin;  other 
dimensions  joules,  cm,  gm,  and  sec.  NC  and  NR  are 
each  one  greater  than  the  number  of  columns  and 
rows.  A  is  the  radius  and  D  the  thickness,  AT 


the  time  step  and  TREF  the  initial  temperature. 

FL  is  the  incident  flux,  and  ABC  the  fraction  ab¬ 


sorbed. 


Arrays,  describing  the  properties  of  the  elements 
are  computed.  Here,  a  homogeneous  material  is  assumed. 
Gaussian  beam  profile  approximated  by  discretiza¬ 
tion.  For  small  SIG  (o)  only  those  cells  entirely 


within  the  beam  receive  flux,  leading  to  significant 


error  if  DR  comparable  to  SIG. 

Beginning  of  main  computational  cycle. 
Temperature  distribution  printed  every  NPRO  time 
steps,  if  desired  (NPUT  ■  0).  Temperatures  not 
printed  for  other  value  of  NPUT. 

Heat  flux  for  all  but  edge  rows  and  columns. 

Heat  flux  for  all  remaining  cells  computed. 

Hoat  balance  performed. 
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